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Abstract—In this paper we present the development of a mathematical 
model for the numerical simulation of one-dimensional compressible 
transient flows in natural gas distribution networks, based on the method of 
characteristics. This model was developed to validate transient phenomena 
in distribution networks of piped natural gas, which is why it was decided 
to build an experimental network with carbon steel pipe, 140 meters in 
length and two inches in diameter, in order to simulate situations involving 
localized ruptures. To achieve the safety conditions and calibration of the 
instruments, the experiments were initially conducted using compressed 
air, with subsequent use of natural gas. The study included several cases of 
leakage in the experimental network, which provided evidence that the 
results obtained show good agreement with the experimental values, thus 
justifying the use of the model for real cases. 





I. INTRODUCTION 


Mathematics plays an important role in several areas of 
study — for example, economics, engineering in general, 
and biological sciences — because through this knowledge 
it is possible to describe the behavior of some specific 
systems or phenomena in mathematical terms. Most 
mathematical formulations for these phenomena lead to 
rates of change of two or more independent variables, 
translating these formulations into partial differential 
equations. The mathematical modeling of these 
formulations may result in distinct approaches, namely: 
experimental, analytical, and computational, which can be 
evaluated jointly or individually. The experimental 
approach requires a physical model that can facilitate 
studies involving the analysis of the direct or indirect 
measurement of the determinant parameters of the 
evaluated problem, which may, in certain circumstances, 
be impractical due to the costs and time involved. For the 
analytical approach, in most cases an adequate solution 
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cannot be obtained because the mathematical techniques 
available are not always suitable for determining such 
solutions. And finally, the computational approach has 
proven to be an important tool because certain 
simplifications allow one to obtain a consistent 
computational model that can be solved using numerical 
methods. 


This work presents a computational mathematical 
formulation that enables the evaluation of transient events 
in distribution networks of piped natural gas, by applying 
the numerical technique based on the method of 
characteristics, which ensures accurate results with 
consistent computational times. To validate the said 
algorithm, we chose to build a corresponding experimental 
model, through which it became possible to simulate 
leakage conditions located in a straight section of piping. 
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I. PARTIAL DIFFERENTIAL EQUATIONS AND 
THEIR CLASSIFICATION 


In engineering, many conventional problems can be 
described by a partial differential equation (PDE) which 
can be classified as elliptical, parabolic, or hyperbolic, 
depending on the category into which the physical 
phenomenon falls (Lax, 2006). Thus, considering a general 
PDE in the following form: 


a> FD T +d 


0x? Oxdy 
.fdtg=0 (1) 





ap Op p p Oh 
c— +d—+e— +e 
dy? Ox Oy 


in which a, b, c, d, e, f, and g may be functions of the 
independent variables x and y and of the dependent 
variable [Iwhich is defined within a region R of the plane 














“xy”, can consider that: 


A = b? — 4ac (2) 


and the PDEs are classified as three distinct types, namely: 











if [1<0, the equation is said to be elliptic; 














if 1=0, the equation is said to be parabolic; and 

















if ||>0, the equation is said to be hyperbolic. 


In terms of problems involving compressible flows, there 
are two distinct evaluation approaches: stationary flows 
and transient flows. In the general stationary problems, the 
differential equations involved are elliptical, similar to 
Laplace’s equation. Transient problems include the 
temporal variation of the magnitudes involved, which is 
why they are represented by hyperbolic or parabolic 
differential equations, corresponding to the wave equation 
and the diffusion equation, respectively. The solution for 
these systems can be achieved either by numerical 
methods or by analytical methods — the method of 
characteristics is both a numerical and analytical method, 
according to Rodrigues (2010). This method is classically 
used in the solution of hyperbolic equations, being a good 
alternative for solving this type of problem, which is why 
it was chosen for the development of this present work. 


Since the PDEs that model physical systems usually have 
many solutions, it is necessary to impose auxiliary 
conditions which may characterize a function that 
represents the solution to the physical problem. Such 
auxiliary conditions correspond to the boundary conditions 
and the initial conditions of the problem. 
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II. MATHEMATICAL FORMULATION FOR THE 


FLOW OF THE GAS 


The transportation of the natural gas along the pipeline 
can be analyzed as follows (Lurie, 2008): the flow is 
compressible and transient; the flow is continuous, causing 
the whole cross-section of the pipe to be filled; the flow 
can be considered to be one-dimensional (i.e., all the 
parameters involved depend only on the x coordinate 
measured along the axis of the pipe and the time t); the 
cross-sectional area of the pipe can vary along the length; 
and the piping can be considered to be indestructible, with 
the interaction between the fluid and the pipe due to 
vibration problems being negligible. 


The mathematical models for fluids and gas flows 
along the pipes are based on physics principles (mechanics 
and thermodynamics) of the continuum and are obtained 
from the following fundamental principles: conservation of 
mass, conservation of momentum, conservation of energy, 
and a corresponding equation of state. 


3.1 — Conservation of mass 


The conservation of mass, or continuity equation, in the 
context of mass flow along a pipe, corresponds to the 
condition in which the mass of the fluid considered can 
neither be created nor destroyed. This condition affects the 
fact that the mass accumulation rate within the control 
volume (CV), outlined in Figure 1, is equal to the net mass 
flow in the control surface; that is: 


Mass flow rate in the | Mass efflux rate inthe |. Mass accumulation rate 
control volume control volume ™) within the control volume 


Fig. 1: Control volume with variable area in a straight 
pipe section (Almeida et al., 2013) 











A At+dA 
p | CV +d P 
——p y De = re v+dV =p 
flow T | T+dT flow 
P i ptdp 
Mathematically, one can write: 
0p 0p OV 10A _ (3) 
ae a a an 0 


given: A = sectional area of the duct; 
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v = velocity of the fluid; 
p = pressure of the fluid; 
T = temperature of the fluid; 


p = fluid density. 


3.2 - Conservation of momentum 


The principle of the conservation of momentum 
corresponds to the application of Newton’s second law of 
motion to a fluid element. Thus, the net force acting, in the 
x direction, on the gas within the control volume, 
corresponds to the algebraic sum of the individual forces 
present in this system in relation to the same reference 
volume (Figure 2); that is: 


M oe Net rate of momentum 
Sum of forces — omentum variation rate = | due to fluid flow across 
within the control volume the boundary 


Fig. 2: Forces acting on the control surface (Almeida et 


al., 2013) 
W 
pone e ee ! 
i —_ | TwÂw | 
i 
e~ me ee — 
( pA). (pA), Hdx 
1 Pal T ww : 
pdf i” 
= dx =| 


For the present case, due to the symmetry of the control 
volume it becomes possible to neglect the radial force 
components, thereby featuring only the presence of forces 
due to the fluid pressure and shear forces (friction). 
Electromagnetic and electrostatic forces are negligible and 
are consequently ignored. Mathematically, this results in: 


ðv 10p Ov 


v— + |v] 


2fov 
ae pon” Ox 7 


D 


0 (4) 





given: Avy = lateral area of the duct; 











w = Shear stress; 





D = diameter of the local duct; 
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fp = friction factor. 


Designating the last component of this expression as Fa 
(Darcy friction force), one can also write: 


OV ðv 1əp (5) 
2 —+-—+F = 
oe “on poe k : 


3.3 Conservation of energy 


The equation for conservation of energy is derived 
from the application of the first Law of Thermodynamics 
which, when formulated for an open system in terms of 
rates, corresponds to: 


E (6) 
— — A 
AL Q-—-W+Ah 


that is, the rate of change of the energy within the control 
volume is equal to the rate of heat transfer into the same 


(Q), minus the work production rate via the control surface 


(W) and, furthermore, added to the net flow of the 
stagnation enthalpy (Uh). 














Since the term corresponding to the work production rate 
can be considered to be zero due to the walls of the ducts 
being rigid, one can obtain the conservation of energy 
equation (explained in terms of the entropy) in the final 
form: 


Os Os kRe- (7) 
stva + Fav) 


given: s= entropy; 
k = ratio of specific heats; 
R = the gas constant; 


c = local speed of sound in the fluid; 


q= heat transfer rate per unit mass. 


3.4 Differential equations in matrix form 


By specific simplifications, with the corresponding 
application of the equation of state for the ideal gas case, 
the group of equations (3), (5), and (7) enables us to write 
the following system of differential equations: 
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oe maun a P a J LR a 
at at at ax. Pax ax” 
7 : pvc? dA 
= p(k-Nq@+vh)- 7G 
naan aa nan CE. an re 
at at at pdx "ax ax” 
= —F, 
fa aed a a ar 
at at at ax ax ax” 
kR . 
= zz (q + VF) 
(08) 

or in the matrix form: 

ðU aU _ 09 

ar A a 


where: 


v cp 0 
1 
B=|—- v 0 
p 
0 0 v 
pvc? dA 


p(k — 1)(q + vF,) — 
f= —Fa 


kR . 
T (q + vF,) 


A dx 


The matrix of coefficients B exhibits real eigenvalues, 
which are different to each other, thus characterizing a 
hyperbolic equation set (Godunov, 1994). Among the 
several methods that can be used to solve hyperbolic 
problems, finite difference techniques based on the Lax- 
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Wendroff scheme and the method of characteristics are 
probably the most used. According to Ames (1965), Lax- 
Wendrof finite difference methods are adequate for 
problems that have a well stated solution within the 
calculation domain. Kruglov et al. (1988), on the other 
hand, pointed out that in the case of flows that may exhibit 
abrupt variations in their properties, the method of 
characteristics can be used to obtain accurate solutions. 


The characteristic curves represent the natural 
coordinate system for solving hyperbolic problems. In this 
coordinate system, the equations in partial derivatives 
defined throughout the entire domain can be expressed as 
ordinary differential equations defined along the 
characteristic curves. These curves are obtained by 
solving the equation set, det( B — dx/dt I) = 0, where I 
is the identity matrix. This equation set expresses the 
necessary condition for the existence of singularities in the 
solution of the problem described by Equation (9) 
(Velasquez et al., 1995). Therefore, the equations of the 
characteristic curves are: 


dx 5 
g c 


dx 
~s 
dx 
a vV 

The first two characteristic curves in the above 
equations are called Mach lines, while the third one is the 
trajectory of the fluid particles (path line). The necessary 
condition for Equation (9) to have a solution along the 
characteristic curves is expressed by the socalled 
compatibility equations: 


Mach line compatibility equations: 


k-1 

(dc) mach IR za) mach 

k-l1e 
= a Ep (4S) macn 
n k—1f vcdA 
2 A dx 

T.a 

— 1) za + vF,) + F,| dt 


+ (k 


(10) 


Path line compatibility equation: 
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kR . 
(ds) natn = (q + vFa)dt 
j (11) 


It should be noted that the compatibility equations for 
the Mach lines in Equation (10) are expressed in terms of 
the differentials of v, c, and s. However, these equations 
can be expressed in terms of only two of these differentials 


























by introducing the Riemann variables || and ||, defined 
according to: 








By doing so, the following compatibility equations 
result: 


c 
(dA), =—(dcy), 
CA 
" k—1f vcdA 
2 A dx 
+(k-1)=(q +R) 
C a 


-F dt 
É (12) 
C 
(dP), = — (dc,)p 
CA 
E k—1;7 vcdA iik 
2 Ai 
Tos 
— 1)—(q + vF,) + F,| dt 
(13) 
— 1c, e 
(dc4) path — > g (q + vF,)dt 
j (14) 


In Equations (12), (13), and (14), a variable called entropy 
level (ca) was introduced instead of entropy s. It is related 
to the former by the following definition (Benson et al., 
1964): 





k—1 
c4 = exp ( s) 
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IV. THE COMPUTATIONAL ROUTINE 


The computational model uses two numerical grids — 
one of them is Eulerian and the other is Lagrangian. The 
Eulerian grid was built by dividing the pipe length into (n - 
1) equal parts, thus identifying nnodes, with two of them at 
the corresponding duct ends. Furthermore, mpoints were 
chosen along the duct, with two of them at the duct ends, 
thus defining mfluid particles whose positions define the 
Lagrangian grid. 


From initial conditions, the fluid velocity and 
thermodynamic properties of the fluid are known along the 
pipe, and starting from these data at (x, t), the 
compatibility equations should be integrated along the 
corresponding characteristic curves, thus giving the 
solution at (x, t +t). The time [lt is chosen in order to 


























satisfy the Courant-Fredrich-Lewy stability criterion 
(Courant et al., 1928), which expresses the necessary 
condition so that the largest displacement of a perturbation 
wave does not exceed the distance between neighboring 











nodes | 1x. 





V. THE EXPERIMENTAL SETUP 


Experimental data were obtained from a pilot network 
approximately 140meters long, which was built of carbon- 
steel piping of 2 inches nominal diameter. In this setup, 
hypothetical transient events were generated by simulating 
localized leaks and using both compressed air and natural 
gas. The experimental apparatus included: five pressure 
sensors installed along the extension of the pipeline in 
order to provide records of the pressure variations 
occurring in the course of each test; a system for regulating 
pressure at the entry of the experimental network to 
maintain the pressure of the system within the limits 
previously established; and, for the test with compressed 
air, a reservoir (pressure vessel) installed at the beginning 
of the network. Figure 3 illustrates an isometric schematic 
of this apparatus, identifying its main components, while 
Figures 4 and 5 depicts specific details such as the vessel 
used during the tests, and the main pipeline. 


It is worth mentioning that, in order to check the 
repeatability of the measurements, each test was performed 
at least three times. Although the opening of the valve was 
performed manually and without the use of any control 
device to guarantees the repeatability of this process, all 
the tests sought to make this opening as quickly as 
possible. By doing it this way, a delay of about 0.5 s was 
observed from the moment of activation to the moment 
that the maximum opening of the valve was reached. 
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Fig. 4: The pressure vessel 





Fig. 5: The main pipeline 
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VI. RESULTS AND DISCUSSION 


Prior to the comparison of numerical and experimental 
results, a study was conducted to verify the results 
obtained after a given number of tests performed in the 
field under localized leakage conditions. Such tests were 
initially performed using compressed air as the test fluid 
and then the gas itself derived from the existing local 
distribution network. 


In Figure 6 it is observed that at the beginning the 
velocity is zero throughout the whole duct and, after 
opening the valve located at x = 137.8 m, which simulates 
the rupture that causes the leakage of the air, this velocity 
increases rapidly at the end of the duct, reaching the sonic 
condition at the throat section of the valve (Figure 8). In 
turn, Figure 7 shows that a fall in pressure at the end of the 
valve occurs simultaneously with the increase in velocity, 
thus generating a depression wave that propagates in the 
direction of the tank and reaches it at approximately t = 
0.45 s. As the depression wave sweeps through the length 
of the duct, the local pressure decreases and causes the 
mass of fluid that lies ahead to begin flowing. The 
depression wave generated at the end of the valve is 
reflected at the other end as a compression wave which, 
upon passing through a point, provokes a reduction in the 
rate at which the local pressure decreases due to the gas 
leakage. This compression wave reaches the end of the 
valve at approximately t = 0.8 s and returns to being 
reflected as a depression wave. This process of depression 
and compression wave propagation causes gas leakages to 
occur as jets, the intensity of which decreases over time. 


a 


Fig. 6: Velocity of the fluid during emptying of the vessel 
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Fig. 7: Pressure of the fluid during emptying of the vessel 
and pipeline 


Discharge valve 
Concentric reducer 





0 1 2 3 4 


Time, [s] 


Fig. 8: Mach number at throat sections of the discharge 
valve and concentric reducer 


In Figures 9-10, the experimental data are compared 
with the results of the simulation for the two cases studied. 
As can be seen in these graphs, the calculated values show 
good agreement with the experimental data, despite the 
final ones displaying noise in the data acquisition system. 
It should be noted that, despite not using any control 
device during testing that would guarantee the 
repeatability of the process of opening the discharge valve, 
in the computational model the flow area in the valve was 
considered to vary according to a sinusoidal function and 
the opening is completed in 0.5 s. It is expected that such 
an approach would be one of the sources of discrepancies 
results and the 


observed between the numerical 


experimental data. 
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Pressure, [har] 





Fig. 9: Numerical and experimental pressure data to a 
generic test 
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Fig. 10: Numerical and experimental pressure data to a 
generic test 


VII. CONCLUDING REMARKS 


A specific mathematical model to describe transient 
events in natural gas distribution networks, based on the 
method of characteristics, was presented in this work. 
Additionally, tests were performed in a pilot pipeline 
system in order to compare numerical and experimental 
thus 
encouraging the use of the mentioned method for solving 


data. The comparison showed good agreement, 


problems of this nature. 
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